library(plotly)

uplan<-function(x,alpha, omega, beta, phi, endow, d, t){-((endow-(beta*phi*(endow-x) + (x*d)/(x+beta*phi*(endow-x))))+(alpha/omega)*((x*d)/(x+beta*phi*(endow-x))-x)-(1-alpha)*t*x)}
uplan2<-function(x,alpha,omega){ uplan(x,alpha,omega,phi=0.4, beta=0.8, d=100, t=0.5, endow=100)}
alpha <- seq(0.001, 1, 0.002)
omega <- seq(0.001, 1, 0.002)

# constrained opt
optx2<-  matrix(nrow=length(alpha),ncol = length(omega))
for(i in 1:length(alpha)){for(j in 1:length(omega)){optx2[i,j] <-constrOptim(theta=c(0.1),f=uplan2,NULL,method = "Nelder-Mead",ci=c(0),ui=matrix(c(1)),alpha=alpha[i],omega=omega[j])$par}}
xyz2<-plot_ly(y=alpha,x=omega,z=optx2, type="surface",colorbar = list(title = "Repatriations",tickmode="array", tickvals=c(0.5,10,34), ticktext=c("Zero","Low","High")))%>%
  layout(scene = list(
    xaxis=list(title='Omega'),
    yaxis=list(title='Alpha'),
    zaxis=list(title='Repatriations',tickmode="array", tickvals=c(0.5,34), ticktext=c("Zero","High"))))

xyz2


#t and beta

uplan<-function(x,alpha, omega, beta, phi, endow, d, at){-((endow-(beta*phi*(endow-x) + (x*d)/(x+beta*phi*(endow-x))))+(alpha/omega)*((x*d)/(x+beta*phi*(endow-x))-x)-(1-alpha)*at*x)}
uplan2<-function(x,beta,at){ uplan(x,beta,at,omega=0.4, phi=0.8, d=100, alpha=0.7, endow=100)}
at <- seq(0.001, 1, 0.002)
beta <- seq(0.001, 1, 0.002)

# constrained opt
optx2<-  matrix(nrow=length(beta),ncol = length(at))
for(i in 1:length(beta)){for(j in 1:length(at)){optx2[i,j] <-constrOptim(theta=c(0.6),f=uplan2,NULL,method = "Nelder-Mead",ci=c(0),ui=matrix(c(1)),beta=beta[i],at=at[j])$par}}
xyz2<-plot_ly(y=beta,x=at,z=optx2,type="surface",colorbar = list(title = "Repatriations",tickmode="array", tickvals=c(0.5,4,13), ticktext=c("Zero","Low","High")))%>%
  layout(scene = list(
    xaxis=list(title='t'),
    yaxis=list(title='Beta'),
    zaxis=list(title='Repatriations',tickmode="array", tickvals=c(0.3,13), ticktext=c("Zero","High"))))

xyz2


# boxplots

#grid with values for phi, beta and tau
valuesbox <- c(0.025, 0.05, 0.10,0.125,0.15,0.20,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,0.975)
alphabox  <- c(0.5, 0.5025,0.505,0.51,0.525,0.55,0.575,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.975)

parambox <- matrix(nrow=3375,ncol = 3)
parambox[,1]<-rep(valuesbox,each=225) #eg phi
parambox[,2]<-rep(valuesbox,each=15)  #eg beta
parambox[,3]<-rep(valuesbox,each=1)   #eg tau

#sigma a
uplan<-function(x,alpha, omega, beta, phi, endow, d, at){-((endow-(beta*phi*(endow-x) + (x*d)/(x+beta*phi*(endow-x))))+(alpha/omega)*((x*d)/(x+beta*phi*(endow-x))-x)-(1-alpha)*at*x)}
uplan2<-function(x,alpha,phi, beta, at){ uplan(x,alpha,omega=0.5,phi, beta, d=100, at, endow=100)}
# constrained opt
boxplot1 <-  matrix(nrow=3375,ncol = length(alphabox))
for(colbox in 1:length(alphabox)){for(param in 1:nrow(boxplot1)){boxplot1[param,colbox] <- 
  constrOptim(theta=c(0.1),f=uplan2,NULL,method = "Nelder-Mead",ci=c(0),ui=matrix(c(1)),alpha=alphabox[colbox],phi=parambox[param,1],beta=parambox[param,2],at=parambox[param,3])$par}}


boxplot.matrix(boxplot1,xlab="Alpha",ylab="Optimal repatriations",names=alphabox,range=0
               ,border=c("blue"))

#sigma b   beta*(phi*endow- max(0,x-(1-phi)*endow))
uplan<-function(x,alpha, omega, beta, phi, endow, d, at){-((endow-(beta*(phi*endow- max(0,x-(1-phi)*endow)) + (x*d)/(x+beta*(phi*endow- max(0,x-(1-phi)*endow)))))+(alpha/omega)*((x*d)/(x+beta*(phi*endow- max(0,x-(1-phi)*endow)))-x)-(1-alpha)*at*x)}
uplan2<-function(x,alpha,phi, beta, at){ uplan(x,alpha,omega=0.5,phi, beta, d=100, at, endow=100)}
# constrained opt
boxplot2 <-  matrix(nrow=3375,ncol = length(alphabox))
for(colbox in 1:length(alphabox)){for(param in 1:nrow(boxplot2)){boxplot2[param,colbox] <- 
  constrOptim(theta=c(0.1),f=uplan2,NULL,method = "Nelder-Mead",ci=c(0),ui=matrix(c(1)),alpha=alphabox[colbox],phi=parambox[param,1],beta=parambox[param,2],at=parambox[param,3])$par}}


boxplot.matrix(boxplot2,xlab="Alpha",ylab="Optimal repatriations, sigma b",names=alphabox,range=0
               ,border=c("blue"))


#sigma c   beta*max(0,phi*endow-x)
uplan<-function(x,alpha, omega, beta, phi, endow, d, at){-((endow-(beta*max(0,phi*endow-x) + (x*d)/(x+beta*max(0,phi*endow-x))))+(alpha/omega)*((x*d)/(x+beta*max(0,phi*endow-x))-x)-(1-alpha)*at*x)}
uplan2<-function(x,alpha,phi, beta, at){ uplan(x,alpha,omega=0.5,phi, beta, d=100, at, endow=100)}
# constrained opt
boxplot3 <-  matrix(nrow=3375,ncol = length(alphabox))
for(colbox in 1:length(alphabox)){for(param in 1:nrow(boxplot3)){boxplot3[param,colbox] <- 
  constrOptim(theta=c(0.1),f=uplan2,NULL,method = "Nelder-Mead",ci=c(0),ui=matrix(c(1)),alpha=alphabox[colbox],phi=parambox[param,1],beta=parambox[param,2],at=parambox[param,3])$par}}


boxplot.matrix(boxplot3,xlab="Alpha",ylab="Optimal repatriations, sigma c",names=alphabox,range=0
               ,border=c("blue"))


# positive solutions as function of phi

#grid with values for phi, beta and tau
valuesbox <- c(0.025, 0.05, 0.10,0.125,0.15,0.20,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,0.975)

parambox <- matrix(nrow=3375,ncol = 3)
parambox[,1]<-rep(valuesbox,each=225) #eg tau
parambox[,2]<-rep(valuesbox,each=15)  #eg beta
parambox[,3]<-rep(valuesbox,each=1)   #eg phi

uplan<-function(x,alpha, omega, beta, phi, endow, d, at){-((endow-(beta*phi*(endow-x) + (x*d)/(x+beta*phi*(endow-x))))+(alpha/omega)*((x*d)/(x+beta*phi*(endow-x))-x)-(1-alpha)*at*x)}
uplan3<-function(x,alpha,omega,phi, beta, at){ uplan(x,alpha,omega,phi, beta, d=100, at, endow=100)}
# constrained opt
possolphi <-  array(numeric(0),dim=(c(length(valuesbox),length(valuesbox),3375)))
for(param in 1:nrow(parambox)){for(alphai in 1:length(valuesbox)){for(omegai in 1:length(valuesbox)){possolphi[alphai,omegai,param] <- 
  constrOptim(theta=c(0.1),f=uplan3,NULL,method = "Nelder-Mead",ci=c(0),ui=matrix(c(1)),alpha=valuesbox[alphai],omega=valuesbox[omegai],phi=parambox[param,3],beta=parambox[param,2],at=parambox[param,1])$par}}}

possolphi2 <-  numeric(3375)
for(j in 1:3375)(possolphi2[j]<-100*2*sum(possolphi[,,j]>0.0001)/225)


dataplot_phi<- data.frame(x=rev(valuesbox), 
                          y1=rev(possolphi2[3361:3375]),
                          y2=rev(possolphi2[1981:1995]),
                          y3=rev(possolphi2[1681:1695]),
                          y4=rev(possolphi2[1231:1245]),
                          y5=rev(possolphi2[2431:2445]),
                          y6=rev(possolphi2[2371:2385]),
                          y7=rev(possolphi2[961:975]),
                          y8=rev(possolphi2[3016:3030]),
                          y9=rev(possolphi2[541:555])
)

plot(dataplot_phi$x,dataplot_phi$y1,type='l',xlab="phi",ylab="% of combinations of alpha and omega with X>0")
lines(dataplot_phi$x,dataplot_phi$y2)
lines(dataplot_phi$x,dataplot_phi$y3)
lines(dataplot_phi$x,dataplot_phi$y4)
lines(dataplot_phi$x,dataplot_phi$y5)
lines(dataplot_phi$x,dataplot_phi$y6)
lines(dataplot_phi$x,dataplot_phi$y7)
lines(dataplot_phi$x,dataplot_phi$y8)
lines(dataplot_phi$x,dataplot_phi$y9)


#endogenous capital controls graphs
alpha <- seq(0.001, 1, 0.01)
omega <- seq(0.001, 1, 0.01)
optendphi2<-  matrix(nrow=length(alpha),ncol = length(omega))
optendx2<-  matrix(nrow=length(alpha),ncol = length(omega))
for(i in 1:length(alpha)){for(j in 1:length(omega)){
  
  endcc<-function(x,alpha,omega, beta, d, t, cmax) {
    return (-(cmax*sqrt(x[1])-(1-alpha)*t*x[2]+(alpha/omega)*((d*x[2]/(x[2]+beta*x[1]*(cmax*sqrt(x[1])-x[2])))-x[2])
              -(beta*x[1]*(cmax*sqrt(x[1])-x[2])+(d*x[2]/(x[2]+beta*x[1]*(cmax*sqrt(x[1])-x[2])))))) }
  eval_ecc2 <- function(x)
  {
    return ( endcc(x,alpha=alpha[i],omega=omega[j], beta=0.8, d=100, t=0.9, cmax=100))  
  }
  # Inequality constraints
  eval_g_ineq <- function(x)
  {
    return (- x[2])
  }
  # Inequality constraints
  eval_g_ineq <- function (x) {
    constr <- c(- x[2]
                
    )
    return (constr)
  }
  # Lower and upper bounds
  lb <- c(0,-0.5)
  ub <- c(1,100)
  #initial values
  x0 <- c(0.2,0.1)
  
  # Set optimization options.
  local_opts <- list( "algorithm" = "NLOPT_LD_MMA", "xtol_rel" = 1.0e-09 )
  opts <- list( "algorithm"= "NLOPT_GN_ISRES",
                "xtol_rel"= 1.0e-09,
                "maxeval"= 100000,
                "local_opts" = local_opts,
                "print_level" = 0 )
  
  res <- nloptr ( x0 = x0,
                  eval_f = eval_ecc2,
                  lb = lb,
                  ub = ub,
                  eval_g_ineq = eval_g_ineq,
                  opts = opts
  )
  
  optendphi2[i,j] <-res$solution[1]
  optendx2[i,j] <-res$solution[2]
}
}

xyz_x2<-plot_ly(y=alpha,x=omega,z=optendx2,type="surface",showscale=FALSE )%>%
  layout(scene = list(
    yaxis=list(title='Alpha'),
    xaxis=list(title='Omega'),
    zaxis=list(title='Repatriations')))
xyz_phi2<-plot_ly(y=alpha,x=omega,z=optendphi2,type="surface",showscale=FALSE )%>%
  layout(scene = list(
    yaxis=list(title='Alpha'),
    xaxis=list(title='Omega'),
    zaxis=list(title='Degree of capital controls')))

xyz_x2
xyz_phi2
